/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include "blockopt.h"

#include "em.h"
#include "fastonebigheader.h"
#include "qfunc.h"

#include <string.h>

float
blockopt (uint64_t                               t,
          float                                  eta,
          float                                  rho,
          uint64_t*                              last_t,
          unsigned int                           n_workers,
          float*                                 dalpha,
          unsigned int                           n_items,
          float*                                 alpha,
          unsigned int                           n_labels,
          float*                                 logbeta,
          float*                                 gamma,
          float*                                 pz,
          const float*                           logpriorz,
          int                                    clamp,
          const Rating*                          ratings,
          unsigned int                           n_ratings,
          const ParameterizedNormalDistribution* prioralpha,
          const Distribution*                    priorlogbeta,
          bool                                   symmetric,
          bool                                   test_only)
{
  /* In the typical crowdsourcing setup a worker can only work on 
   * an item once, so don't try to guard against duplicate workers.
   */

  /* Step 1: apply alpha and gamma prior sparsely 
   * NB: this is just a (cheap!) approximation
   */

  /* NB: Limit rho->1 exists, but I'm lazy. */

  if (! test_only)
    {
      float powt = fastpow (t, 1.0f - rho);
      float sigma = prioralpha->sq_stddev; 

      for (unsigned int i = 0; i < n_ratings; ++i)
        {
          uint64_t s = last_t[ratings[i].worker];
          float pows = fastpow (s, 1.0f - rho);
          float decay = 
            fastexp (eta * (pows - powt) / (n_items * (1.0f - rho) * sigma));

          float* alphaptr = alpha + ratings[i].worker * n_labels * n_labels;

          for (unsigned int j = 0; j < n_labels; ++j)
            {
              for (unsigned int k = 0; k < n_labels; ++k)
                {
                  if (k == j)
                    {
                      /* do nothing */
                    }
                  else if (symmetric && k < j)
                    {
                      float dalphatmp = 
                        (1.0 - decay) * 
                        (gamma[k * n_labels + j] - alphaptr[k * n_labels + j]);

                      alphaptr[k * n_labels + j] += dalphatmp;
                      gamma[k * n_labels + j] += dalphatmp / (1.0 + n_workers);
                    }
                  else
                    {
                      float dalphatmp = 
                        (1.0 - decay) * 
                        (gamma[j * n_labels + k] - alphaptr[j * n_labels + k]);

                      alphaptr[j * n_labels + k] += dalphatmp;
                      gamma[j * n_labels + k] += dalphatmp / (1.0 + n_workers);
                    }
                }
            }

          last_t[ratings[i].worker] = t;
        }
    }

  /* Step 2: maximize per-item objective wrt logbeta and pz */

  float q = em (alpha, n_labels, logbeta, pz, logpriorz, clamp, 
                ratings, n_ratings, priorlogbeta, symmetric, 100, 1e-3);

  if (! test_only)
    {
      /* Step 3: gradient update of alpha and gamma */

      for (unsigned int i = 0; i < n_ratings; ++i)
        {
          memset (dalpha + ratings[i].worker * n_labels * n_labels, 
                  0, 
                  n_labels * n_labels * sizeof (dalpha[0]));
        }

      dqfunc 
        (dalpha, alpha, n_labels, *logbeta, pz, ratings, n_ratings, symmetric);

      float etat = eta * fastpow (t, -rho);

      for (unsigned int i = 0; i < n_ratings; ++i)
        {
          float* alphaptr = alpha + ratings[i].worker * n_labels * n_labels;
          float* dalphaptr = dalpha + ratings[i].worker * n_labels * n_labels;

          /* this is where duplicate workers would cause multiple updates */

          for (unsigned int j = 0; j < n_labels * n_labels; ++j)
            {
              alphaptr[j] += etat * dalphaptr[j];
              gamma[j] += (etat / (1.0 + n_workers)) * dalphaptr[j];
            }
        }

      /* priorgamma is assumed Gaussian with unit variance
       * and has already been incorporated into step 1
       */
    }

  return q;
}

void
finalize (uint64_t                               t,
          float                                  eta,
          float                                  rho,
          uint64_t*                              last_t,
          unsigned int                           n_workers,
          unsigned int                           n_items,
          float*                                 alpha,
          unsigned int                           n_labels,
          float*                                 gamma,
          const ParameterizedNormalDistribution* prioralpha,
          bool                                   symmetric)
{
  /* Step 1: apply deferred alpha and gamma prior
   * NB: this is just a (cheap!) approximation
   */

  /* NB: Limit rho->1 exists, but I'm lazy. */

  float powt = fastpow (t, 1.0f - rho);
  float sigma = prioralpha->sq_stddev; 

  for (unsigned int w = 0; w < n_workers; ++w)
    {
      uint64_t s = last_t[w];
      float pows = fastpow (s, 1.0f - rho);
      float decay = 
        fastexp (eta * (pows - powt) / (n_items * (1.0f - rho) * sigma));

      for (unsigned int j = 0; j < n_labels; ++j)
        {
          for (unsigned int k = 0; k < n_labels; ++k)
            {
              if (k == j)
                {
                  /* do nothing */
                }
              else if (symmetric && k < j)
                {
                  float alphasave = 
                     alpha[w * n_labels * n_labels +
                           k * n_labels + 
                           j];

                  float dalphatmp = 
                    (1.0 - decay) * (gamma[k * n_labels + j] - alphasave);

                  alpha[w * n_labels * n_labels + 
                        k * n_labels + 
                        j] += dalphatmp;

                  gamma[k * n_labels + j] += dalphatmp / (1.0 + n_workers);
                }
              else
                {
                  float alphasave = 
                     alpha[w * n_labels * n_labels +
                           j * n_labels + 
                           k];

                  float dalphatmp = 
                    (1.0 - decay) * (gamma[j * n_labels + k] - alphasave);

                  alpha[w * n_labels * n_labels + 
                        j * n_labels + 
                        k] += dalphatmp;

                  gamma[j * n_labels + k] += dalphatmp / (1.0 + n_workers);
                }
            }
        }

      last_t[w] = t;
    }
}
